Nearly-one-dimensional self-attractive Bose-Einstein condensates in optical lattices 
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Within the framework of the mean-field description, we investigate atomic Bose-Einstein conden- 
sates (BECs) , with attraction between atoms, under the action of strong transverse confinement and 
periodic (optical-lattice, OL) axial potential. Using a combination of the variational approximation 
(VA), one- dimensional (ID) nonpolynomial Schrodinger equation (NPSE), and direct numerical so- 
lutions of the underlying 3D Gross-Pitaevskii equation (GPE) , we show that the ground state of the 
condensate is a soliton belonging to the semi-infinite bandgap of the periodic potential. The soliton 
may be confined to a single cell of the lattice, or extend to several cells, depending on the effective 
self-attraction strength, g (which is proportional to the number of atoms bound in the soliton), 
and depth of the potential, Vb, the increase of Vb leading to strong compression of the soliton. We 
demonstrate that the OL is an effective tool to control the soliton's shape. It is found that, due to 
the 3D character of the underlying setting, the ground-state soliton collapses at a critical value of the 
strength, g — g c , which gradually decreases with the increase of the depth of the periodic potential, 
Vb; under typical experimental conditions, the corresponding maximum number of 7 Li atoms in the 
soliton, -/V max , ranges between 8,000 and 4,000. Examples of stable multi-peaked solitons are also 
found in the first finite bandgap of the lattice spectrum. The respective critical value g c again slowly 
decreases with the increase of Vb, corresponding to 7V max ~ 5,000. 

PACS numbers: 03.75.Lm; 05.45.Yv; 03.75.Hh 



I. INTRODUCTION 

It has been firmly established in experiments with 
ultracold vapors of alkali metals that Bose-Einstein 
condensates (BECs) with weak attractive interactions 
between atoms (characterized by negative scattering 
length) can form stable matter-wave solitons in nearly 
one-dimensional (ID) "cigar-shaped" traps, which tightly 
confine the condensate in two transverse directions, but 
leave it almost free along the longitudinal axis. This set- 
ting has made it possible to create stable bright solitons 
0, IH and trains of such solitons [l[ in the 7 Li conden- 
sate, where the interaction between atoms may be made 
weakly attractive by means of the Feshbach-resonance 
technique. In the 85 Rb condensate trapped under sim- 
ilar conditions, stronger attraction between atoms leads 
to the creation of nearly 3D solitons in a post-collapse 
state p|. 

This experimentally relevant setting is described by 
effectively ID equations which were derived, by means 
of various approximations, from the full 3D Gross- 
Pitaevskii equation (GPE) [l]-[l5j]. In some cases, the 
deviation of the effective equation from the straightfor- 
ward one-dimensional GPE with the cubic nonlinearity 
is adequately accounted for by an extra self-attractive 
quintic term, which may essentially affect the solitons 
0) B Q • The derivation of the ID equation starts with 
adopting an ansatz factorizing the 3D wave function into 
the product of a transverse one (which amounts to the 
ground state of the 2D harmonic oscillator) and slowly 
varying axial (one-dimensional) wave function. The sub- 
stitution of this ansatz in the variational approximation 



(VA; a review of the method can be found in Ref. 
leads, without resorting to additional approximations, to 
a nonpolynomial Schrodinger equation (NPSE) for the 
longitudinal wave function [H, Q . This equation was suc- 
cessfully used in various physical contexts [n]> El- I n 
particular, the above-mentioned simplified equation in- 
cluding cubic and quintic terms can then be obtained 
by an expansion of the NPSE for the case of a relatively 
weak nonlinearity A generalization of the NPSE for a 
two-component BEC was recently developed in Ref. [l3| , 
in the form of a coupled system of NPSEs. 

The objective of this work is to analyze the dynamics 
of quasi- ID solitons in the cigar-shaped trap equipped 
with a periodic potential, which can be easily created as 
an optical lattice (OL), by means of two coherent coun- 
terpropagating laser beams illuminating the condensate 
(for a recent review of the BEC dynamics in periodic 
potentials, see Ref. (HI)- Solitons in the ID nonlin- 
ear Schrodinger equation (NLS) with the ordinary self- 
focusing cubic nonlinearity and periodic (lattice) poten- 
tial had been studied some time ago [l7|. It was shown 
that a stable soliton can be trapped by a local potential 
minimum of the OL. On the other hand, stable quasi- 
1D solitons in cigar-shaped traps exist up to a certain 
threshold, beyond which they suffer collapse. The occur- 
rence of the collapse in the effective ID equation reflects 
thepresence of the collapse in the underlying 3D setting 
U& [fa l2l • A new feature, which the present paper aims to 
report and explore, is the influence of the periodic poten- 
tial on the collapse threshold for solitons in the quasi- ID 
trap. We demonstrate that, for both single-peak solitons 
found in the semi-infinite gap of the periodic potential, 
and multi-peaked solitons found in finite bandgaps, the 
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collapse threshold gradually goes down with the increase 
of the potential strength. This feature may be explained 
by the fact that a sufficiently strong axial OL potential 
collects almost all the condensate in a limited space of a 
single lattice cell, thus facilitating the onset of the col- 
lapse. The predicted effect should be amenable to obser- 
vation by means of standard experimental techniques. 

The paper is organized as follows. In Section II, we 
apply the VA directly to the description of solitons in 
the underlying 3D GPE. In this way, the collapse thresh- 
old for the quasi-lD solitons is predicted in a semi- 
analytical form [in the same approximation, the stability 
of the solitons is estimated by means of the Vakhitov- 
Kolokolov (VK) criterion]. In Section III, we derive the 
one-dimensional NPSE in the presence of the periodic 
potential, and find soliton families as numerical solutions 
of the latter equation. Comparison with direct numerical 
solutions of the underlying 3D GPE demonstrates that 
the NPSE provides for very high accuracy in the predic- 
tion of both the shape of the solitons belonging to the 
semi-infinite gap and their collapse threshold; the accu- 
racy provided by the VA is lower, but nevertheless rea- 
sonable too. In Section IV, we briefly consider examples 
of solitons belonging to the first finite bandgap. They are 
found in a numerical form from the corresponding NPSE. 
In the semi-infinite gap, the solitons may be narrow, oc- 
cupying, essentially, a single site of the lattice potential, 
or broad, extending to several sites, but they always fea- 
ture a single tall peak. The soliton in the finite bandgap 
has a very different shape, with many peaks. The paper 
is concluded by Section V. 



the transverse harmonic length, aj_ = \/h/ (mwi) (m is 
the atomic mass), and the depth of the potential, Vo, is 
taken in units of Huj±, therefore uj± does not appear in 
the equations. For 7 Li atoms in the transverse trap with 
•jj±_ = 2t: x 1 KHz, one has a± ~ 1 ^tm; then, = 1 (the 
value used in exampled displayed below) and A = 1.5 
/im correspond to 9 ~ 30°. In the same case, Vo = 1 is 
tantamount, in physical units, to ~ 1.5E rcc , where the 
recoil energy for the OL with period d = A/ (2sin(#/2)) 
is E rec — (nh/d) /m. 

To predict solitons in an approximate analytical form, 
we use the 3D Gaussian ansatz, 
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where a and r\ are, respectively, the transverse width and 
axial length of the localized pattern. Inserting this ansatz 
into Eqs. (fTJ) and @, we obtain 
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Aiming to predict the ground state in the framework 
of the above approximation, we look for values of a and 
t] that minimize energy E [as given by Eq. ([5])]. using 
equations dE/da = dE/drj = 0. This way, we derive 
coupled equations, 



II. VARIATIONAL APPROACH 
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The energy-per-atom (E) and chemical potential (fj,) 
of the self-attractive BEC described by the mean-field 
stationary wave function, ip(r), in the presence of the 
strong transverse harmonic confinement with frequency 
acting in the plane of (x,y), are 
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Here, the OL potential acting along axis z is 
U(z) = -Vbcos(2fc L z) , 
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with Ul — (2n/X) sin (9/2), where A is the wavelength 
of two laser beams with angle 9 between them that cre- 
ate the OL. We assume normalization J \tp(r)\ 2 dr = 1, 
then g = 2\a s \N/a± is the adimensional strength of the 
self-attraction, with negative scattering length of atomic 
collisions a s , and the number of atoms in the condensate, 
N. Lengths in Eqs. fTJ) - ([3]) are measured in units of 
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which can be solved numerically [18j. Obviously, these 
solutions yield a ground state, i.e., a minimum of energy, 
only if the curvature of the energy dependence, E(r],a), 
dispositive, 
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As concerns the dynamical stability of solitons against 
small perturbations, it may be, first of all, estimated by 
means of the VK criterion [l9( , according to which a nec- 
essary stability condition is d[ij dg < (in the present no- 
tation), if the soliton family is described by dependence 
fi(g). Actually, the VK criterion may also be sufficient for 
the stability of solitons in the GPE with self-attraction 
and lattice potential [20| . 

In Fig. 1, we plot solutions for a and rj of Eqs. (0) and 
([5]), together with the corresponding energy-per-particle 
E [calculated as per Eqs. (J5J|], as functions of interaction 
strength g, for several fixed values of potential depth Vq 
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FIG. 2: Continuation of Fig. [T] chemical potential n of the 
quasi-lD soliton versus g, as predicted by the variational ap- 
proximation for the same three values of Vo, and kh — 1. 



FIG. 1: (Color online). Axial length r\, transverse width a, 
energy-per-particle E of the quasi-lD soliton versus the self- 
attraction strength, g = 2\a 3 \N/a±, as predicted by the vari- 
ational approximation based on ansatz Q and Eqs. Q, ©■ 
The dependences are displayed at fixed values of depth Vo 
of periodic potential @, with k^ = 1. On the right side, all 
curves terminate at a critical point, g = g c (Vo), beyond which 
the collapse occurs (the upper curve in the top panel, which 
is cut by the panels's frame, actually extends up to g = 0, 
like its counterpart in the bottom panel. 



and wavenumber hh = 1. Additionally, Fig. 2 displays 
respective dependences fi(g) of the soliton's chemical po- 
tential, found from Eq. ©. The figures show that the 
soliton in the self-attractive BEC is predicted to exist 
up to a critical strength, g c , which depends on Vo and 
k,L, At g > g c , the 3D collapse of the nearly- ID soli- 
ton occurs, as suggested by results obtained previously 
in models without the OL potential [J, [H, Q . Note that 
the VA predicts g c = 1.55 for Vo — , which is somewhat 
higher than the critical value, g c — 1.33, obtained from a 
numerical solution of the full GPE with Vq = in three 
dimensions @, [24[ • The characteristic value of the non- 
linearity strength in Figs. 1 and 2, g = 1, corresponds, 
for the above-mentioned transverse size, a± = 1 /im, and 
scattering length a s = —0.1 nm, which can be attained in 
7 Li by means of the Feshbach resonance [H, , to solitons 
built of N ~5, 000 atoms. 

As g drops to zero, the axial size of the soliton, 
77, diverges, while transverse width a approaches 1 (in 
physical units, it becomes equal to the above-mentioned 
harmonic-oscillator length, a±). On the other hand, as 
g approaches g Cl both rj and a remain finite and smaller 
than 1. 



At small Vo, the figure shows only a small distortion 
of the curves, and a small reduction of g c , in comparison 
with the previously studied case of Vo — 0. A qualitative 
change [which is best pronounced in Fig. 2, in terms 
of the ji(g) dependences] is observed at Vo > 0.15 (i.e., 
for > 1,000 atoms, according to the above estimate), 
when there appear two different stable branches. For 
V = 0.2 and Vo = 0.4, the lower branches of the 77(17) and 
cr(g) dependences exists only in a finite interval, which 
we denote as 



9m < 9 <9c 
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while the upper branches extend up to g — 0, i.e., they 
exist in interval < g < gM, with gM < g c - Physically, 
the lower branches (with smaller values of axial length 
77) correspond to an attractive BEC which is localized, 
essentially, within a single cell of the OL; accordingly, 
we call the corresponding solution a single-site soliton. 
The upper branches of 77(g) and a(g) correspond to a 
weakly localized solution, which occupies several lattice 
sites; therefore, we call it a multi-site soliton. The bot- 
tom panel of Fig. 1 demonstrates that, for Vo > 0.15, 
the multi-site soliton provides for smaller energy, i.e., 
it represents the ground state, at small g. At larger g, 
the single-site soliton becomes the ground state, while its 
multi-site counterpart is a metastable state. In particu- 
lar, the above estimates demonstrate that, for Vo = 0.4 
(which is tantamount to ~ 0.6-Brec), the switch from the 
multi-site ground state to the single-site one occurs when 
the number of atoms (in the 7 Li condensate) attains val- 
ues ~ 2,500. 

Figure 2 suggests that the families of soliton solutions 
are dynamically stable, as they always meet the VK cri- 
terion, d/i/dg < 0. In the case when two solutions exist, 
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FIG. 3: The critical value of the effective self-attraction 
strength, g c , above which the soliton collapses, versus wave 
number fcz, of periodic potential Q. The results are obtained 
from the numerical solution of variational equations and 



III. NONPOLYNOMIAL SCHRODINGER 
EQUATION 

A more accurate analysis of the present setting may 
be performed using the one-dimensional nonpolynomial 
Schrodinger equation (NPSE) [5(. As mentioned in In- 
troduction, the derivation of the NPSE starts with the 
ansatz factorizing the 3D wave function into the prod- 
uct of the transverse 2D one, which describes the ground 
state of the harmonic oscillator, and a slowly varying ID 
axial function, f(z) (which may be complex), 



1 



exp 



{x 2 + V 2 ) 
2a{z) 2 



(10) 



Obviously, it is a more general ansatz than the above one, 
based on Eq. (|4]). Inserting expression (p~0|) into Eq. (Q}, 
one finds [upon neglecting spatial derivatives of &(z)] the 
following one-dimensional energy functional 
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i.e., at Vo > 0.15, both of them are VK-stable, i.e., the 
ground state and its metastable counterpart alike. 

For Vq > 0.45, the numeral solution of Eqs. {7} and 
([5]) demonstrates that the lower threshold g m of existence 
interval © for the single-site soliton vanishes, but this is 
as an artifact of the VA, which is actually known in other 
contexts too (2(j. It is explained by inaccuracy of the 
Gaussian ansatz in the limit of weak nonlinearity, i.e., for 
widely spread small-amplitude solitons. In this situation, 
one may, in principle, apply a more sophisticated ansatz, 
combining the Gaussian and periodic functions, such as 
cos(2fciz); however, the generalized ansatz results in a 
cumbersome algebra [2l| , therefore we do not follow this 
way here. 

It is interesting to predict the collapse critical strength, 
<? c , as a function of parameters Vo and /cl of the axial pe- 
riodic potential. In Fig. 3, we display dependence g c (ki,) 
predicted by the VA at four fixed values of Vo. For given 
strength Vo of the potential, there exists wavenumber k c 
at which g c attains its minimum, i.e., the collapse has the 
lowest threshold. Further, Fig. 3 shows that this mini- 
mum of g c decreases with the increase of Vo , which may 
be understood as mentioned above: the strong poten- 
tial tends to squeeze the entire condensate into a single 
call of the lattice, which facilitates the onset of collapse. 
On the other hand, at large values of /cl, g c is asymp- 
totically constant, as the interaction of the condensate 
with the short-period OL becomes exponentially weak, 
see Eq. ([5]), hence it produces little effect on the collapse 
threshold. 
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Further, imposing an extra normalization condition, 
dz\f(z)\ 2 = l, 
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on the ID wave function and adding a term with the cor- 
responding Lagrange multiplier, /i (which will again be 
the chemical potential), to energy functional (fTTjl . and, 
finally, minimizing the energy, one arrives at the follow- 
ing equations for real functions f(z) and a(z) [where the 
expression ([3]) for the potential is substituted]: 
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Equation (flU)) is the stationary version of the NPSE || 
with the periodic potential, while /i is fixed by normaliza- 
tion condition (fT2|) . Note that only the ID wave function, 
f(z), appears in the NPSE, while the transverse width is 
locally expressed in terms of f(z), as per Eq. (fT4|) . Equa- 
tion (fT"3"|) reduces to the familiar ID cubic NLS equation 
for gf 2 <C 1 [then, Eq. yields a w 1], Only in 

the latter case, the system may be considered as truly 
one-dimensional. 

In previous studies which used the NPSE [l| [(| [ll|, 
GJ, [H, solitons were not considered in the presence of 
the periodic potential. In this work, we solved the full 
one-dimensional NPSE, which includes the time deriva- 
tive and the periodic potential, numerically, using the 
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FIG. 4: The axial density profile, \f(z)\ 2 , of the soliton in 
periodic potential with fcz, = 1 and four different values 
of Vq. The self-attraction strength is fixed at g = 0.5. 



finite-difference Crank-Nicholson method in imaginary 
time [22j. In this way, we have obtained the soliton pro- 
files shown in Fig. 4 for different values of the poten- 
tial depth Vo and fixed self-attraction strength, g = 0.5 
(recall it typically corresponds to ~ 2, 500 atoms of 7 Li 
bound in the soliton). As seen from the figure, the in- 
crease of the OL strength from zero to ~ 2 recoil energies 
may compress the soliton, in the axial direction, by a fac- 
tor > 4. Generally, the use of the OL offers an efficient 
means to control the soliton shape (the OL may also be 
used as versatile tool to manipulate matter-wave solitons 
dynamically [HI). 

For Vo = (dashed line in Fig. 1), the shape of the 
soliton has a unique maximum, and the density profile 
may be well fitted to f(z) = (y/g/2) sech(gz/2), which 
is rigorously valid in the above-mentioned limit corre- 
sponding to gf 2 <C 1, provided that z varies in infinite 
limits. On the other hand, if z belongs to a finite interval, 
—L/2 < z < L/2, with periodic boundary conditions, 
f(z + L) = f(z), the NPSE with V = yields a spatially 
uniform profile of the ground state, f[z) = l/VX, for 
sufficiently weak nonlinearity, < g < it 2 /L; t he g round 
state develops a spatial structure at g > it 2 / L (l2Tl2a|. 

For nonzero but small Vb (dot-dot-dashed and solid 
curves in Fig. 3), the soliton profile features several lo- 
cal maxima and minima due to the effect of the periodic 
potential [2a] . Thus, under such conditions, the Bosc 
condensate self-traps into a multi-peaked soliton, which 
occupies several cells of the periodic potential (the ex- 
istence of stable three-dimensional multi-peaked solitons 
with a similar shape in the periodic axial potential, but 
without any transverse confinement, was demonstrated 
in Ref. [271 ]: in that case, the transverse self-trapping 
and stability of the solitons was provided for by the 
"Feshbach-resonance-management" technique, i.e., peri- 
odic alternation of the sign of nonlinear coefficient g). 
Following the increase of Vo, the soliton compresses in 
the axial direction, and (for instance, at Vq = 1.2, see 
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FIG. 5: Axial length of the ground-state bright soliton, 
(z 2 ) 1//2 , as a function of self-attraction strength g, for Vo = 

0. 4 and fez, = 1 in Eq. @- Displayed are results provided by 
the variational approximation based on the Gaussian ansatz, 

1. e., Eqs. J7| and (JSj , and by the nonpolynomial Schrodinger 
equation (NPSE), Eq. (TT5)) . 



the dot-dashed curve in Fig. 4) the secondary maxima 
become very small. In this case, the condensate actually 
self-traps into a single-peak soliton, which occupies only 
one cell of the OL ©. 

Contrary to the sharp transition predicted by the 
VA in the previous section, the NPSE shows a smooth 
crossover from the multi-peak soliton to the single-peak 
one. In Fig. 5, we plot the axial length, (z 2 ), of the 
ground-state bright soliton as a function of strength g, 
for Vo = 0.4 (which is tantamount to the OL strength 
~ 0.6 recoil energy, as shown above), and compare the 
prediction of the VA with results produced by the NPSE. 
It is clearly seen that, while the VA predicts a jump at 
g ~ 0.5, the NPSE does not show it. In fact, this jump 
is a consequence of the fact that ansatz (Q}, which as- 
sumes the simple Gaussian waveform for the axial wave 
function, provides for a crude fit to multi-peaked states. 
We also note that the numerical solution of the NPSE 
does not reveal the bistability predicted by the VA, as 
the numerical algorithm seeks for the ground-state solu- 
tion for given g. For typical values of physical parameters 
mentioned above, Fig. 5 implies that, as the number of 
atoms increases from 1, 000 to 5, 000, the soliton shrinks 
in the longitudinal direction from 10 to 1 /xm. 

The VA developed in the previous section predicts the 
collapse of the soliton at the critical value of the self- 
attraction strength, g = g c , which depends on parame- 
ters Vb and fez, of periodic potential ([3]). It is relevant to 
compare this prediction with results following from nu- 
merical solution of NPSE (fT3")l . In Table 1, we present 
the critical value, g c , found from the NPSE for different 
values of Vb and = 1. The results are in qualita- 
tive agreement with their variational counterparts dis- 
played in Fig. 3: the increase of the potential depth, 
Vq, leads to gradual reduction of g c . Note that the value 
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of g c = 1.33 for Vq = is virtually the same as the 
above-mentioned numerically exact one, g c = 1.34, which 
was found from the numerical solution of the full three- 
dimensional GPE (with V = 0) dEi]. According to the 
above estimates, the drop of g c from 1.33 to 0.85 corre- 
sponds to the decrease of the maximum number of atoms 
from N max ~ 8, 000 to 7V max ~ 4, 000. For the sake of 
completeness, the table also shows the axial length and 
transverse width, obtained from Eqs. (fl3)) and (fl"4f at 
the collapse point, g — g c . It is seen that the soliton 
shrinks in both directions with the increase of Vq, but 
its length and widths remain finite up to the collapse 
point. Typical values of the physical parameters referred 
to above imply that the data presented in Table 1 pre- 
dict roughly constant density of the condensate at the 
collapse threshold, ~ 3 x 10 14 cm -3 . 
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Table 1. The critical value of the self-attraction 
strength, g c , and the corresponding values of the axial 
length, \J (z 2 ), and minimal transverse width, <x(0), of the 
soliton in the periodic potential, V(z) = — Vq cos [2k lz), 
with fci = 1, for different values of Vq, as found from 
numerical solution of the ID nonpolynomial Schrodinger 
equation, Eq. (TT3]). 

An important issue is the comparison of the results 
yielded by the NPSE with those found from direct nu- 
merical solution of the full 3D GPE in the presence of 
the axial periodic potential (in previous works dealing 
with the NPSE |g, El, EH, El, such a comparison was 
not presented for the model including the OL). It is also 
interesting to compare the results with those which can 
be obtained from the ordinary ID cubic GPE [i.e., the 
ID equation with the cubic nonlinearity and the same 
periodic potential, which can be obtained from Eq. (|13p 
by expanding in small g and keeping only terms up to 
the first order in g] , since the latter equation is frequently 
used as a model of the BEC in the quasi-lD traps. Figure 
6 shows that the density profiles generated by the NPSE 
are always very close to (practically, coincide with) the 
ones obtained from the 3D GPE. On the other hand, the 
profiles generated by the ID cubic GPE are different, and 
the difference gets more pronounced with the increase of 
self- attraction strength g. This observation is not sur- 
prising because the nonlinearity in the NPSE essentially 
deviates from the cubic term if g\f(z)\ 2 is large enough. 
Note also that the cubic GPE in one dimension cannot 
predict any collapse. 
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FIG. 6: The axial density profile, p(z), of the soliton in po- 
tential @, with fc_L = 1 and Vo = 0.2. Comparison between 
results provided by the different equations: NPSE, 3D GPE, 
and ID GPE (the latter two with the cubic nonlinearity) is 
presented. In the case of the 3D equation, Eq. (TTJ), the ax- 
ial density is defined as p(z) = J J \ip(r)\ 2 dxdy, while in the 
other cases it is simply |/(z)| 2 . 



IV. SOLITONS IN THE FINITE BANDGAP 

The above analysis was dealing with solitons whose 
chemical potential belongs to the semi-infinite bandgap 
in the linear spectrum of Eq. (fl~3"J) . On the other hand, it 
is known that the cubic self-attractive nonlinearity may 
also give rise to solitons located in higher-order (finite) 
bandgaps [28| . This section addresses such solitons in the 
present model. It is relevant to stress that up to now, 
higher-bandgap solitons, in the case of self-attraction, 
have been considered only in strictly one-dimensional set- 
tings. 

Here, we report soliton solutions found by means of the 
NPSE, Eq. (|13|) . in the first finite bandgap correspond- 
ing to ki, — 1. For this purpose, a self-consistent numeri- 
cal method was used, with periodic boundary conditions. 
We employed a spatial grid of 1025 points, covering the 
interval of —50.26 < z < 50.26, which corresponds to 32 
periods of the external potential. To verify the correct- 
ness of the numerical scheme, we have checked that the 
lowest-energy state in the semi- infinite bandgap (i.e., the 
ground state of the system), produced by this method, 
is identical to that found above by the integration of the 
NPSE in imaginary time, i.e., the soliton shown in Fig. 
4. 

As a typical example, in Fig. 7 we plot the density 
profile, |/(z)| 2 , of the 33th state (it is number 33 in the 
full set of the states generated by the numerical scheme) 
for g = 0.5 and three different values of the potential 
depth, Vq. In the linear approximation, this state lies at 
the bottom of the second Bloch band. With the increase 
of nonlinearity strength g, the energy of the 33th state 
lowers, and, in doing so, it enters the first finite bandgap 
from above. At large values of g (and small values of 
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FIG. 7: The axial-density profile, |/(z)| , of the lowest state in 
the first finite bandgap of periodic potential ((3} with kz, = 1. 
The self-attraction strength is g = 0.5. The top, middle, and 
bottom panels correspond, respectively, to Vb = 0.2, Vb = 0.6, 
and Vb = 1. 



strates gradual decrease of gi X \ The latter dependence 
is qualitatively the same as reported above for solitons in 
the semi- infinite bandgap, cf. Table 1. 

The axial size of the soliton at the collapse threshold 
(.9 = Sc 1 **) is also included in Table 1. As for the trans- 
verse width, a(z), in all cases it takes value <r(0) = 1 at 
the center (z = 0), unlike the solitons in the semi- infinite 
gap, cf. Table 1. This feature is explained by the fact 
that the amplitude of the gap soliton vanishes at z = 
(see Fig. 7), hence the width of the confined state at this 
point is the same as in the linear equation, i.e., a = 1 
(according to the normalization adopted above). 

In Table 2, we start with Vq — 4 because for smaller 
Vb the considered state plunges into the semi-infinite 
bandgap at g close to gc . According to the above es- 
timates, the largest number of Li atoms possible in the 
soliton created in the first bandgap is N max ~ 5, 000. 

For the sake of completeness, in Fig. 8 we present 
the first 41 eigenvalues e,-, as found from the stationary 
NPSE equation, 



Vb), it passes the entire first bandgap from its top to the 
bottom, then crosses the band separating this bandgap 
from the semi-infinite gap, and eventually sinks into the 
latter one. Figure 7 shows that, for Vq = 0.2, this state 
is still fully delocalized (being similar to a Bloch wave), 
while for Vq = 0.6 it becomes localized, with a mean- 
square width much smaller than the length of the periodic 
box, featuring many local maxima and minima (zeros). 
We call this solution an excited multi-site soliton. Its 
multi-peaked structure is typical to gap solitons [HI, [H, 
[29j . At Vb = 1, Fig. 7 shows that the excited soliton 
compresses itself into a narrower state. Simulations of the 
NPSE in real-time demonstrate that this soliton family 
is dynamically stable. 



1 d 2 



n/1-s]/;(z)I 2 



Vb 


o (1) 




4 


1.16 


0.70 


5 


1.10 


0.63 


6 


1.03 


0.60 


7 


1.00 


0.57 



Table 2. The critical value of the self- attraction 
strength, g^\ of the lowest-energy state in the first fi- 
nite bandgap of periodic potential ([3]), with hi, = 1 and 
increasing values of its depth, Vq. At g > g ^, c ollapse 
takes place. The soliton's axial length, y/ (z 2 ), corre- 
sponding to to g = <7c , is tabulated too. Note that 
<r(0) is always equal to 1 because the soliton axial profile 
|/(z)| 2 has a node at z = 0. 

For the lowest-energy state in the first finite bandgap, 
there also exists a critical strength, g^\ of the self- 
attraction strength, above which the collapse occurs. In 
Table 2, we present these critical values corresponding to 
the increasing depth of the potential, Vb, which demon- 



(15) 



and plotted versus interaction strength g. The respective 
eigenfunctions, fj(z), may be both delocalized (Bloch- 
like) and localized (soliton- like). At g = 0, the first 32 
eigenvalues generated by our numerical scheme belong to 
the first band, while the other nine fall into the second 
band. In compliance with the above discussion, Fig. 8 
shows that the lowest eigenvalue and the 33th one split 
off from the first and second continuous bands and move 
down (up to the onset of the collapse) with the increase of 
<7, thus giving rise to localized states, in the semi-infinite 
and first finite gaps, respectively. It is noteworthy that 
the second and third eigenvalues, which originally belong 
to the first band, also split off from it at larger values of 
g (at largest values of g displayed in Fig. 8, the fourth 
eigenvalues demonstrates the same trend). We have veri- 
fied that the corresponding nonlinear eigenstates become 
localized, as one may expect. Note that similar results 
have been found in Ref. [3C| in a numerical solution of 
the ordinary cubic GPE in one dimension. 



V. CONCLUSIONS 

In this work, we have considered soliton states in self- 
attractive BECs loaded into a cigar-shaped trap, which is 
equipped with a periodic OL (optical-lattice) axial poten- 
tial. Using the 3D variational approximation (VA), and 
the effective NPSE (nonpolynomial Schrodinger equa- 
tion) in one dimension, we have demonstrated that the 
ground state of the condensate is a bright soliton whose 
chemical potential falls within the semi-infinite bandgap 
created by the periodic potential. With the increase of 
the self-attraction strength, g, and/or the OL strength, 
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FIG. 8: Eigenvalues ej found from Eq. (fTo)) with Vo = 1 and 
fcz, = 1, as functions of self-attraction strength g. 

Vo, the ground-state soliton changes its shape from broad 
(multi-site) to a narrow (single-site) one. In particular, 
the soliton composed of ~ 2, 500 7 Li atoms trapped in 
the transverse harmonic potential with oj± = 2ir x 1 KHz 
gets compressed by a factor of 4, as Vq increases from zero 
to two recoil energies; generally, the OL offers a conve- 
nient tool to control the shape of the soliton. The soliton 
solutions produced by the VA are stable against small 
perturbations, according to the VK (Vakhitov-Kolokolov 
criterion). 

Due to the 3D nature of the trapped self-attractive 
condensate, the solitons exist up to a critical value, g c , of 
g, beyond which the collapse takes place. A new feature 
reported in this work is the gradual decrease of g c with 
the increase of Vo , a qualitative explanation to which was 
given. For typical values of experimentally relevant pa- 
rameters, the results translate into an estimate for the 
critical number of atoms at which the collapse occurs, 



which ranges from 8, 000 to 4, 000, while the critical den- 
sity of the collapsing condensate is, roughly, constant, 
~ 3 x 10 14 cm~ 3 . Comparison with numerical solutions of 
the underlying 3D Gross-Pitaevskii equation shows that 
the NPSE very accurately predicts both the shape of the 
ground-state solitons and dependence g c (Vo); the accu- 
racy of the variational approximation is somewhat lower, 
but also reasonable. 

Stable multi-peaked solitons were also found in the first 
finite bandgap of the OL potential. The respective criti- 
cal nonlincarity strength again slowly decreases with the 
increase of Vo, the corresponding largest number of atoms 
in the soliton being estimated as ~ 5,000. It was stud- 
ied in detail how the localized state splits off from the 
Bloch band and gives rise to the soliton in the first finite 
bandgap. 

Finally, it is relevant to mention that an interest- 
ing possibility, which deserves further consideration, is 
to induce an effective axial potential not directly, but 
rather through periodic modulation of the transverse 
tight-binding frequency (l4| . Another promising exten- 
sion of the present analysis may be to the case of self- 
repulsive BEC loaded in the combination of the cigar- 
shaped trap and axial periodic potential. The study of 
the latter setting may shed new light on the theoretical 
description of quasi-lD gap solitons [lj|[3l|]. These issues 
will be considered elsewhere. 
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